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An analytical model of unsteady heat transfer in a one-dimensional harmonic crystal is presented. 
A nonlocal temperature is introduced as a generalization of the kinetic temperature. A closed 
equation determining unsteady thermal processes in terms of the nonlocal temperature is derived. 
For an instantaneous heat perturbation a time-reversible equation for the kinetic temperature is 
derived and solved. The resulting constitutive law for the heat flux in the considered system is 
obtained. This law significantly differs from Fourier’s law and it predicts a finite velocity of the heat 
front and independence of the heat flux on the crystal length. The analytical results are confirmed 
by computer simulations. 
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An understanding of heat transfer at microlevel is es¬ 
sential to obtain link between microscopic and macro¬ 
scopic description of solids [1, 2]. As far as macroscopic 
scale level is concerned the Fourier law of heat conduction 
is widely and successfully used to describe heat transfer 
processes. At microscopic level, however, analytical and 
numerical investigations have shown substantial devia¬ 
tions from Fourier’s law [3-5]. These inadequacies can be 
on principle addressed by using special laws of particles 
interactions [6-9] or complex enough structures [10, 11]. 
Recent experimental data however showed that Fourier’s 
law is indeed violated in low-dimensional nanostruc¬ 
tures [12-14]. This motivates interest to the simplest 
lattice models, in particular harmonic one-dimensional 
crystals (chains), where these anomalies are most promi¬ 
nent [15, 16]. Problems of this kind previously have been 
mainly addressed in the context of the steady-state heat 
conduction [3-5]. The present work focuses on unsteady 
conduction regimes [11, 17-19]. 

Here we suggest an approach that allows rigorous 
derivation of macroscopic heat conduction equations and 
corresponding anomalous heat conduction law for har¬ 
monic systems in a one-dimensional, non-quantum case. 
The obtained equations differ substantially from the ear¬ 
lier suggested heat transfer equations [20, 21], however 
they are in excellent agreement with molecular dynam¬ 
ics simulations and previous analytical estimations [18]. 
The approach presented here is developed in context of 
one-dimensional systems, however the same ideas can be 
applied to systems of higher dimensions. 

Thus we consider a one-dimensional crystal, described 
by the following equation of motion: 

Ui = - 2ui + Ui+i), Wo = v^cym) (1) 

where ui is the displacement of the ith particle, m is the 
particle mass, C is the stiffness of the interparticle bond. 
The crystal is infinite: the index i is an arbitrary integer. 
The initial conditions are 

Ui\t=o = ( 2 ) 


where Qi are independent random values with zero expec¬ 
tation and unit variance; cr^(x) is variance of the initial 
velocities, which is a slowly varying function of the spatial 
coordinate x = ia, where a is the lattice constant. These 
initial conditions correspond to an instantaneous temper¬ 
ature perturbation, which can be induced in crystals, for 
example, by an ultrashort laser pulse [22]. The displace¬ 
ments as functions of time Ui = Ui{t) can be found as 
a solution of the Cauchy problem (l)-(2). In addition, 
these functions are linearly dependent on the integration 
constants, which are random due to the random nature 
of the initial conditions (2). 

The first analytical solution of a steady heat conduc¬ 
tion problem for a harmonic chain was obtained in [23] 
using a covariance matrix for coordinates and momenta. 
Then this approach was extended and applied to various 
harmonic systems [3, 10, 15, 16]. Study of the covari¬ 
ance matrix allowed obtaining analytical expressions for 
steady [24] and unsteady [25, 26] temperature profiles. In 
this paper a somewhat similar approach based on analysis 
of covariances for velocities [27-29] is used. We introduce 
a nonlocal temperature 0n{x) as 

ksi-lT 0n{x) m{uiUj), (3) 

where is the Boltzmann constant, n = j — i is the co- 
variance index, x = is the spatial coordinate, angle 
brackets stand for mathematical expectation, {uiUj) is 
the velocity covariance (note that {iii) = {uj) =0). If 
n = 0 then i = j and quantity 0^ coincides with the ki¬ 
netic temperature T: 0o{x) = T{x) = where 

i = xja. According to its definition, the introduced non¬ 
local temperature reflects a nonlocal nature of thermal 
processes in harmonic crystals and can be considered as 
a generalization of the kinetic temperature. 

Let us calculate the second time derivative of On us¬ 
ing the dynamics equations (1) and the following two 
approximations. 

1. The nonlocal temperature 0n{x) is a slowly varying 
function of the spatial coordinate x (on the dis- 
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tances of order of the lattice constant a). This 
allows replacing the finite differences by spatial 
derivatives [30]. The approximation is adequate for 
processes that are sufficiently smooth in space, e. g. 
for spatial temperature profiles in a form of waves 
that are much longer then a. 

2. The virial approximation [1]: time or spatial deriva¬ 
tives of mathematical expectations are small with 
respect to quantities that have non-zero values in 
thermodynamic equilibrium. In particular, this al¬ 
lows us to express covariances of the bond strains 
in terms of the nonlocal temperature. The approx¬ 
imation is adequate for processes that are not too 
far from thermodynamic equilibrium. 

Then after some transformations we obtain a differential- 
difference equation 


of Bessel functions. In particular, Bessel functions were 
successfully applied to solution of shock-wave problems in 
harmonic chains [31, 32]. Similarly, the problem ( 6 ) has 
an analytical solution 0^(1, k) = To(/c) J 2 n(cH), where 
J 2 n are the Bessel functions of the 1st kind [33]. From 
the practical point of view the most interesting case is 
n = 0, which gives Fourier image T{t,k) of the kinetic 
temperature distribution: 

f{t,k)=fo{k)Jo{ckt). (7) 

From (7) it follows that the image T(t^k) satisfies the 
Bessel differential equation 

( 8 ) 

Fourier inversion of ( 8 ) gives a partial differential equa¬ 
tion for the temperature field 


~ — 0, (4) 

where c = ujoa is the speed of sound. This is a closed 
equation describing unsteady thermal processes in the 
crystal in terms of the nonlocal temperature. The pro¬ 
cesses under consideration should be such that the nonlo¬ 
cal temperature is sufficiently smooth in time and space. 
Apart from this limitation any unsteady thermal pro¬ 
cesses satisfy equation (4). After solution (analytical or 
numerical) of equation (4) the kinetic temperature can 
be obtained as T(t,x) = 0n{t,x)\n=o. 

The initial conditions for equation (4) corresponding 
to the original initial conditions ( 2 ) are: 

0n\t=o = To{x)6n, 4 |t =0 = 0 , (5) 

where To{x) = ^^ma‘^{x) is the initial temperature dis¬ 
tribution; = 1 for n = 0 and = 0 for n 7 ^ 0. The 
initial conditions (5) are taken after a fast transition pro¬ 
cess, which results, according to the virial theorem, in a 
double reduction of the initial kinetic temperature [28]. 
Note that in contrast with the random initial value prob¬ 
lem (l)-( 2 ), the initial value problem (4)-(5) is expressed 
in terms of mathematical expectations, and therefore it 
is a deterministic problem. 

Using an integral Fourier transform in the spatial coor¬ 
dinate X the problem (4)-(5) can be solved analytically. 
For the Fourier image 0n{t, k) we obtain 

On = \c^k‘^{0n-l - ‘^On + ^n+l), 

A 

On\t=0 = TQ{k)Sn , On\t=0 = 0 , 

where k is the spatial frequency, To{k) is the Fourier im¬ 
age of the initial temperature distribution To{x). Let us 
note the similarity between (l)-( 2 ) and ( 6 ): initial value 
problem ( 6 ) can be interpreted as a motion of a har¬ 
monic chain having an initial shift of the central particle. 
This kind of problems can be effectively solved in terms 


f + U = c^T". 


(9) 


The corresponding initial conditions follow from (5): 

T\t=o = To{x), t\t=o = 0. ( 10 ) 


Fourier inversion of the representation (7) gives an ana¬ 
lytical solution of the initial value problem (9)-(10): 



To (a; - ct) 

-. (IT. 


( 11 ) 


Similar integral representation was obtained in [34] using 
heat energy density correlation functions. 

Thus, the evolution of the temperature field in a one¬ 
dimensional crystal after an instantaneous thermal per¬ 
turbation is described by partial differential equation (9) 
with initial conditions ( 10 ) or by integral formula ( 11 ). 
According to ( 11 ) the thermal front propagates with the 
sound speed c (in contrast to the thermal conductivity 
based on Fourier’s law where an unphysical instantaneous 
signal propagation is realized). The obtained wave be¬ 
havior of the heat front is similar to predictions of the 
wave theories of heat conduction [20, 21]. However, the 
obtained solution has important differences, which will 
be shown in the text to follow. 

Let us consider the heat flux. For the considered sys¬ 
tem it can be represented [4, 5, 35] as 

q= ^C{{ui-Ui+i){ui+Ui+i)). ( 12 ) 


The heat flux q satisfies the energy balance equation 

pksf = -q\ (13) 


where p = 1 /a is the density (number of particles per unit 
volume). Joint consideration of equations (9) and (13) 
gives the constitutive law for the heat flux 

q+jq = -kspc^T', ( 14 ) 
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which replaces Fourier’s law in the considered system. 
Alternatively, the law (14) can be derived directly, in the 
same way as equation (9) is derived. An integral repre¬ 
sentation for the heat flux follows from (11) and (14): 

Comparison of the obtained heat transfer equation (9) 
with the heat equation based on Fourier’s law and the 
thermal wave equation based on the Maxwell-Cattaneo- 
Vernotte (MCV) law is given in table 1. The latter equa- 



Heat 

equation 

(Fourier) 

Thermal wave 
equation 
(MCV) 

Equation (9) 
(present paper) 

a) 

II 

rjn ^ 1 rjn _ rjn// 

T T 

f -h -T — (At" 
t 

b) 

1 

II 

q-\- -q = --T 

r T 

4 + = -kBpAT' 

c) 


^ cos (kct) 



solution for the temperature and heat flux 


T(t, x) = AoJo{kct) sin kx + 5, 

(17) 

q{t,x) = —kBpcA^Jiikct) coskx, 

where Jq and Ji are the Bessel functions of the 1st kind. 
Previously, an existence of a Bessel function solution for 
this problem was mentioned in [18], and solution similar 
to (17) for the temperature field was obtained in Master- 
degree thesis [38]. 

To justify the assumptions in derivation of the ana¬ 
lytical solution we compare it with results of molecular 
dynamics (MD) simulations. Equations (1) are solved by 
the central differences method, the time step is O.OItq, 
where tq = 27r/ci;o- The initial conditions (16) are set by 
a random number generator, the wavelength A is equal to 
the length of the chain containing 10^ particles. To pro¬ 
vide correspondence with the analytical approach used 
above, 10^ realizations of such chain with an independent 
random initiation are computed. To optimize the com¬ 
putations all chains are joined at end-points to form a 
long chain (10^ particles) with periodic boundary condi¬ 
tions. The results of the computations are compared with 
analytical solution (17) in Fig. 1. The horizontal axis 


TABLE I. a) Heat transfer equation, b) equation connecting 
heat flux and temperature, c) decay law for the sinusoidal 
heat perturbation. Notations: t is time (variable), r is the 
relaxation time (constant), [3 is the thermal diffusivity, k, is 
the thermal conductivity, c is the sound speed, p is the den¬ 
sity, ks is the Boltzmann constant, k is the spatial frequency. 
Approximation (c) for the thermal wave equation is obtained 
for = P/t and large k] approximation for Jo is valid for 
relatively large t. 

tion and equation (9) have similar form and somewhat 
similar behavior (e. g. a finite velocity of the heat front 
propagation). However, there are significant differences. 
The first one is that r, a material constant, is replaced 
in (9) by the physical time t. The second difference is 
time-reversibility of equation (9): the equation is not 
changing when t is replaced by — t, same as the origi¬ 
nal microscopic equation (1). On contrary, both classical 
and wave equations of heat conduction are irreversible. 
The contradiction between time-reversibility of the clas¬ 
sical microscopic equations and irreversibility of the cor¬ 
responding macroscopic continuum equations is one of 
the opened questions of the modern physics [1, 36]. The 
obtained reversible macroscopic equation of heat conduc¬ 
tion may be a step towards solution of this problem. 

We consider now a sinusoidal temperature perturba¬ 
tion [11, 37]: 

Tq{x) = Aosinkx (16) 

where Aq and B are temperature constants, k = 27 rlX 
is the spatial frequency, A is the wavelength of the per¬ 
turbation. Eormulas (11) and (15) give exact analytical 



FIG. 1. Oscillational decay of the thermal perturbation am¬ 
plitude for ID harmonic crystal. Comparison of the analytical 
solution (17) with the MD modeling results (10^ joined chains 
containing 10^ particles each). Dashed lines show the enve¬ 
lope proportional to 1/v^ and also an exponential envelope 
inherent to the thermal wave equation. 

in Fig. 1 represents the dimensionless time t/to, where 
to = A/c; the vertical axis stands for the oscillation am¬ 
plitude A(t), which is computed as the first coefficient of 
a spatial Fourier expansion of the temperature field. Ac¬ 
cording to Fig. 1 there is an excellent agreement between 
the computational results and the analytical curve. 

Due to the Bessel function properties [33], the tem¬ 
perature and heat flux (17) have an oscillational decay, 
where the oscillation amplitude is asymptotically pro¬ 
portional to l/Vi. The same asymptotics has been ob¬ 
tained in [18] for one-dimensional harmonic crystals. In 
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Fig. 1 the envelope proportional to is shown by 

the dashed lines, perfectly bounding both analytical and 
computational graphs. The existing theories of heat con¬ 
duction [20, 21], such as Fourier’s, Maxwell-Cattaneo- 
Vernotte (MCV), dual-phase-lag [39], and spacetime- 
elasticity [19] yield linear differential equations with con¬ 
stant coefficients, and therefore all of them predict an ex¬ 
ponential decay of the sinusoidal perturbation amplitude. 
In table 1 a comparison of the analytically obtained de¬ 
cay law for A{t)/AQ with the results based on some other 
theories is demonstrated, an exponential envelope inher¬ 
ent to the thermal wave model is also shown in Fig. 1. 
Thus, among the mentioned theories only the current one 
gives an analytical solution, which agrees with the MD 
simulations and asymptotic estimations of the oscillation 
decay for harmonic chains [18]. 

Let us consider now a stepwise initial temperature dis¬ 
tribution, modeling heat transfer between a hot and a 
cold body: 

x<0: T(x)=T 2 , x>0: T{x)=Ti, (18) 

where T 2 > Ti. In this case the integral representations 
(11), (15) yield for \x\ < ct an exact analytical solution 

T{t,x) =Ti + ^arccos^, 

^_ (19) 

where AT = T 2 — Ti; for x > ct the original tempera¬ 
ture distribution remains and the heat flux is zero. Ac¬ 
cording to (19) the heat front propagates with the sound 
speed c and the heat flux through cross-section x = 0 
is constant and equal to ^kspcAT. In contrast, use of 
Fourier’s law for the same problem gives the heat flux 
proportional to which is infinite at t = 0 (an un¬ 

physical consequence of Fourier’s law). Thus the heat 
flux ^kBpcAT is provided by the temperature difference 
that is realized on the spatial interval x G [—ct, ct] with 
increasing length of 2ct. Consequently, the heat flux de¬ 
pends on the temperature difference rather than on the 
temperature gradient. This is in qualitative agreement 
with the known phenomenon of thermal superconductiv¬ 
ity: the heat flux through a one-dimensional harmonic 
crystal placed between two thermal reservoirs does not 
depend on the length of the crystal [5, 23]. The same 
value ^kspcAT was obtained in [40] as a steady-state 
limit of the heat flux for large t. 

In Fig. 2 the analytical solution (19) is compared with 
computer simulations for T 2 = 2Ti. The above described 
computation procedure is used. Fig. 2 shows the ini¬ 
tial temperature distribution, the analytical solution, and 
the computation results obtained at t = to/8 using 10^ 
and 10^ particles (to = L/c^ where L is the chain length; 
only half of the chain is shown in the figure). Conver¬ 
gence to the analytical solution with the increase of the 
system size is clearly seen. 



FIG. 2. Heat transfer between hot (left) and cold (right) 
areas of ID harmonic crystal. The analytical solution (19) 
is compared with the computer simulation (MD): 10^ chains 
containing 10^ particles each (cross is an average over 10 par¬ 
ticles); 10^ chains containing 10^ particles (circle is an average 
over 100 particles). 


Fig. 3 shows a part of Fig. 2 corresponding to pos¬ 
itive X. For symmetry reasons this case can be inter¬ 
preted as a problem of a half-space heating: the initial 
temperature for x > 0 is Ti and the boundary condi¬ 
tion at X = 0 is T = (T 2 + Ti)/2 > Ti. The advantage 



FIG. 3. Heat propagation for different ID models: a) heat 
equation, b) thermal wave equation, c) wave equation, 
d) equation (9), e) computer simulation for ID harmonic crys¬ 
tal (10^ chains containing 10^ particles each). 

of this formulation is that the constant boundary tem¬ 
perature is maintained without any thermostat. This is 
important since the heat transfer can substantially de¬ 
pend on the thermostat properties [41, 42]. Solutions of 
the considered problem using four different continuum 
equations are compared in Fig. 3 with the simulation re¬ 
sults. Parameters are chosen in such a way that the total 
heat quantity transferred through the cross-section x = 0 
(area under each curve) is equal for all models and the 
heat front (when it exists) propagates with the sound 
speed c. According to Fig. 3 the computation results al¬ 
most coincide with the analytical solution of equation (9) 
and significantly differ from the solutions based on the 
other theories of thermal conduction. Indeed, the clas¬ 
sical heat equation predicts no heat front, the thermal 
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wave equation gives a stepwise front, while the real heat 
front is described by a smooth curve having a vertical 
tangent at x = ct. Note that the thermal wave (MCV) 
equation behaves as wave equation at small times and 
as heat equation at large times [43]. However, according 
to the analytical solution (19) and the presented com¬ 
puter simulations, the heat transfer in a one-dimensional 
harmonic crystal is self-similar, i. e. T = T(^). 

Thus, a nonlocal temperature (a generalization of the 
kinetic temperature) is introduced in this work to ob¬ 
tain closed description of thermal transfer in a one¬ 
dimensional harmonic crystal. Finally this yields to a 
partial differential equation (9) for the kinetic tempera¬ 
ture, which can be referred to as the time-reversible ther¬ 
mal wave equation. The resulting macroscopic constitu¬ 
tive law (14) (an alternative of Fourier’s law for the con¬ 
sidered system) predicts a finite velocity of the heat front 
and independence of the heat flux on the crystal length. 
The analytical findings are in excellent agreement with 
the molecular dynamics simulations and previous ana¬ 
lytical estimations. The obtained results are relevant to 
aspects of nanotechnology that involve heat transfer pro¬ 
cesses in high purity nanostructures [12, 13, 44]. 
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